<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demgp</title>
  <meta name="keywords" content="demgp">
  <meta name="description" content="DEMGP	Demonstrate simple regression using a Gaussian Process.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">netlab</a> &gt; demgp.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../menu.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\netlab&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>demgp
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>DEMGP	Demonstrate simple regression using a Gaussian Process.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">DEMGP    Demonstrate simple regression using a Gaussian Process.

    Description
    The problem consists of one input variable X and one target variable
    T. The values in X are chosen in two separated clusters and the
    target data is generated by computing SIN(2*PI*X) and adding Gaussian
    noise. Two Gaussian Processes, each with different covariance
    functions are trained by optimising the hyperparameters  using the
    scaled conjugate gradient algorithm.  The final predictions are
    plotted together with 2 standard deviation error bars.

    See also
    <a href="gp.html" class="code" title="function net = gp(nin, covar_fn, prior)">GP</a>, <a href="gperr.html" class="code" title="function [e, edata, eprior] = gperr(net, x, t)">GPERR</a>, <a href="gpfwd.html" class="code" title="function [y, sigsq] = gpfwd(net, x, cninv)">GPFWD</a>, <a href="gpgrad.html" class="code" title="function g = gpgrad(net, x, t)">GPGRAD</a>, <a href="gpinit.html" class="code" title="function net = gpinit(net, tr_in, tr_targets, prior)">GPINIT</a>, <a href="scg.html" class="code" title="function [x, options, flog, pointlog, scalelog] = scg(f, x, options, gradf, varargin)">SCG</a></pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="gp.html" class="code" title="function net = gp(nin, covar_fn, prior)">gp</a>	GP	Create a Gaussian Process.</li><li><a href="gpcovar.html" class="code" title="function [cov, covf] = gpcovar(net, x)">gpcovar</a>	GPCOVAR Calculate the covariance for a Gaussian Process.</li><li><a href="gpfwd.html" class="code" title="function [y, sigsq] = gpfwd(net, x, cninv)">gpfwd</a>	GPFWD	Forward propagation through Gaussian Process.</li><li><a href="gpinit.html" class="code" title="function net = gpinit(net, tr_in, tr_targets, prior)">gpinit</a>	GPINIT	Initialise Gaussian Process model.</li><li><a href="netopt.html" class="code" title="function [net, options, varargout] = netopt(net, options, x, t, alg);">netopt</a>	NETOPT	Optimize the weights in a network model.</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="demnlab.html" class="code" title="function demnlab(action);">demnlab</a>	DEMNLAB A front-end Graphical User Interface to the demos</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">%DEMGP    Demonstrate simple regression using a Gaussian Process.</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%    Description</span>
0004 <span class="comment">%    The problem consists of one input variable X and one target variable</span>
0005 <span class="comment">%    T. The values in X are chosen in two separated clusters and the</span>
0006 <span class="comment">%    target data is generated by computing SIN(2*PI*X) and adding Gaussian</span>
0007 <span class="comment">%    noise. Two Gaussian Processes, each with different covariance</span>
0008 <span class="comment">%    functions are trained by optimising the hyperparameters  using the</span>
0009 <span class="comment">%    scaled conjugate gradient algorithm.  The final predictions are</span>
0010 <span class="comment">%    plotted together with 2 standard deviation error bars.</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%    See also</span>
0013 <span class="comment">%    GP, GPERR, GPFWD, GPGRAD, GPINIT, SCG</span>
0014 <span class="comment">%</span>
0015 
0016 <span class="comment">%    Copyright (c) Ian T Nabney (1996-2001)</span>
0017 
0018 
0019 <span class="comment">% Find out if flops is available (i.e. pre-version 6 Matlab)</span>
0020 v = version;
0021 <span class="keyword">if</span> (str2num(strtok(v, <span class="string">'.'</span>)) &gt;= 6)
0022     flops_works = logical(0);
0023 <span class="keyword">else</span>
0024     flops_works = logical(1);
0025 <span class="keyword">end</span>
0026 
0027 randn(<span class="string">'state'</span>, 42);
0028 x = [0.1 0.15 0.2 0.25  0.65 0.7 0.75 0.8 0.85 0.9]';
0029 ndata = length(x);
0030 t = sin(2*pi*x) + 0.05*randn(ndata, 1);
0031 
0032 xtest = linspace(0, 1, 50)';
0033 
0034 clc
0035 disp(<span class="string">'This demonstration illustrates the use of a Gaussian Process'</span>)
0036 disp(<span class="string">'model for regression problems.  The data is generated from a noisy'</span>)
0037 disp(<span class="string">'sine function.'</span>)
0038 disp(<span class="string">' '</span>)
0039 disp(<span class="string">'Press any key to continue.'</span>)
0040 pause
0041 
0042 flops(0);
0043 <span class="comment">% Initialise the parameters.</span>
0044 net = <a href="gp.html" class="code" title="function net = gp(nin, covar_fn, prior)">gp</a>(1, <span class="string">'sqexp'</span>);
0045 prior.pr_mean = 0;
0046 prior.pr_var = 1;
0047 net = <a href="gpinit.html" class="code" title="function net = gpinit(net, tr_in, tr_targets, prior)">gpinit</a>(net, x, t, prior);
0048 
0049 clc
0050 disp(<span class="string">'The first GP uses the squared exponential covariance function.'</span>)
0051 disp(<span class="string">'The hyperparameters are initialised by sampling from a Gaussian with a'</span>)
0052 disp([<span class="string">'mean of '</span>, num2str(prior.pr_mean), <span class="string">' and variance '</span>, <span class="keyword">...</span>
0053     num2str(prior.pr_var), <span class="string">'.'</span>])
0054 disp(<span class="string">'After initializing the network, we train it using the scaled conjugate'</span>)
0055 disp(<span class="string">'gradients algorithm for 20 cycles.'</span>)
0056 disp(<span class="string">' '</span>)
0057 disp(<span class="string">'Press any key to continue'</span>)
0058 pause
0059 
0060 <span class="comment">% Now train to find the hyperparameters.</span>
0061 options = foptions;
0062 options(1) = 1;    <span class="comment">% Display training error values</span>
0063 options(14) = 20;
0064 flops(0)
0065 [net, options] = <a href="netopt.html" class="code" title="function [net, options, varargout] = netopt(net, options, x, t, alg);">netopt</a>(net, options, x, t, <span class="string">'scg'</span>);
0066 <span class="keyword">if</span> flops_works
0067     sflops = flops;
0068 <span class="keyword">end</span>
0069 
0070 disp(<span class="string">'The second GP uses the rational quadratic covariance function.'</span>)
0071 disp(<span class="string">'The hyperparameters are initialised by sampling from a Gaussian with a'</span>)
0072 disp([<span class="string">'mean of '</span>, num2str(prior.pr_mean), <span class="string">' and variance '</span>, num2str(prior.pr_var)])
0073 disp(<span class="string">'After initializing the network, we train it using the scaled conjugate'</span>)
0074 disp(<span class="string">'gradients algorithm for 20 cycles.'</span>)
0075 disp(<span class="string">' '</span>)
0076 disp(<span class="string">'Press any key to continue'</span>)
0077 pause
0078 flops(0)
0079 net2 = <a href="gp.html" class="code" title="function net = gp(nin, covar_fn, prior)">gp</a>(1, <span class="string">'ratquad'</span>);
0080 net2 = <a href="gpinit.html" class="code" title="function net = gpinit(net, tr_in, tr_targets, prior)">gpinit</a>(net2, x, t, prior);
0081 flops(0)
0082 [net2, options] = <a href="netopt.html" class="code" title="function [net, options, varargout] = netopt(net, options, x, t, alg);">netopt</a>(net2, options, x, t, <span class="string">'scg'</span>);
0083 <span class="keyword">if</span> flops_works
0084     rflops = flops;
0085 <span class="keyword">end</span>
0086 
0087 disp(<span class="string">' '</span>)
0088 disp(<span class="string">'Press any key to continue'</span>)
0089 disp(<span class="string">' '</span>)
0090 pause
0091 clc
0092 
0093 fprintf(1, <span class="string">'For squared exponential covariance function,'</span>);
0094 <span class="keyword">if</span> flops_works    
0095     fprintf(1, <span class="string">'flops = %d'</span>, sflops);
0096 <span class="keyword">end</span>
0097 fprintf(1, <span class="string">'\nfinal hyperparameters:\n'</span>)
0098 format_string = strcat(<span class="string">'  bias:\t\t\t%10.6f\n  noise:\t\t%10.6f\n'</span>, <span class="keyword">...</span>
0099   <span class="string">'  inverse lengthscale:\t%10.6f\n  vertical scale:\t%10.6f\n'</span>);
0100 fprintf(1, format_string, <span class="keyword">...</span>
0101     exp(net.bias), exp(net.noise), exp(net.inweights(1)), exp(net.fpar(1)));
0102 fprintf(1, <span class="string">'\n\nFor rational quadratic covariance function,'</span>);
0103 <span class="keyword">if</span> flops_works
0104     fprintf(1, <span class="string">'flops = %d'</span>, rflops);
0105 <span class="keyword">end</span>
0106 fprintf(1, <span class="string">'\nfinal hyperparameters:\n'</span>)
0107 format_string = [format_string <span class="string">'  cov decay order:\t%10.6f\n'</span>];
0108 fprintf(1, format_string, <span class="keyword">...</span>
0109       exp(net2.bias), exp(net2.noise), exp(net2.inweights(1)), <span class="keyword">...</span>
0110       exp(net2.fpar(1)), exp(net2.fpar(2)));
0111 disp(<span class="string">' '</span>)
0112 disp(<span class="string">'Press any key to continue'</span>)
0113 pause
0114 
0115 disp(<span class="string">' '</span>)
0116 disp(<span class="string">'Now we plot the data, underlying function, model outputs and two'</span>)
0117 disp(<span class="string">'standard deviation error bars on a single graph to compare the results.'</span>)
0118 disp(<span class="string">' '</span>)
0119 disp(<span class="string">'Press any key to continue.'</span>)
0120 pause
0121 cn = <a href="gpcovar.html" class="code" title="function [cov, covf] = gpcovar(net, x)">gpcovar</a>(net, x); 
0122 cninv = inv(cn);
0123 [ytest, sigsq] = <a href="gpfwd.html" class="code" title="function [y, sigsq] = gpfwd(net, x, cninv)">gpfwd</a>(net, xtest, cninv);
0124 sig = sqrt(sigsq);
0125 
0126 fh1 = figure;
0127 hold on
0128 plot(x, t, <span class="string">'ok'</span>);
0129 xlabel(<span class="string">'Input'</span>)
0130 ylabel(<span class="string">'Target'</span>)
0131 fplot(<span class="string">'sin(2*pi*x)'</span>, [0 1], <span class="string">'--m'</span>);
0132 plot(xtest, ytest, <span class="string">'-k'</span>);
0133 plot(xtest, ytest+(2*sig), <span class="string">'-b'</span>, xtest, ytest-(2*sig), <span class="string">'-b'</span>);
0134 axis([0 1 -1.5 1.5]);
0135 title(<span class="string">'Squared exponential covariance function'</span>)
0136 legend(<span class="string">'data'</span>, <span class="string">'function'</span>, <span class="string">'GP'</span>, <span class="string">'error bars'</span>);
0137 hold off
0138 
0139 cninv2 = inv(<a href="gpcovar.html" class="code" title="function [cov, covf] = gpcovar(net, x)">gpcovar</a>(net2, x));
0140 [ytest2, sigsq2] = <a href="gpfwd.html" class="code" title="function [y, sigsq] = gpfwd(net, x, cninv)">gpfwd</a>(net2, xtest, cninv2);
0141 sig2 = sqrt(sigsq2);
0142 fh2 = figure;
0143 hold on
0144 plot(x, t, <span class="string">'ok'</span>);
0145 xlabel(<span class="string">'Input'</span>)
0146 ylabel(<span class="string">'Target'</span>)
0147 fplot(<span class="string">'sin(2*pi*x)'</span>, [0 1], <span class="string">'--m'</span>);
0148 plot(xtest, ytest2, <span class="string">'-k'</span>);
0149 plot(xtest, ytest2+(2*sig2), <span class="string">'-b'</span>, xtest, ytest2-(2*sig2), <span class="string">'-b'</span>);
0150 axis([0 1 -1.5 1.5]);
0151 title(<span class="string">'Rational quadratic covariance function'</span>)
0152 legend(<span class="string">'data'</span>, <span class="string">'function'</span>, <span class="string">'GP'</span>, <span class="string">'error bars'</span>);
0153 hold off
0154 
0155 disp(<span class="string">' '</span>)
0156 disp(<span class="string">'Press any key to end.'</span>)
0157 pause
0158 close(fh1);
0159 close(fh2);
0160 clear all;</pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>